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Abstract 

The impulse response function (IRF) of a localized bolus in cerebral blood flow codes 
important information on the tissue type. It is indirectly accessible both from MR- 
and CT-imaging methods, at least in principle. In practice, however, noise and lim- 
ited signal resolution render standard deconvolution techniques almost useless. Para- 
metric signal descriptions look more promising, and it is the aim of this contribution 
to develop some improvements along this line. 
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1 Introduction 

The detailed knowledge of blood circulation in the brain is beyond doubt an important 
diagnostic tool. Among other aspects, it admits the detection and a partial characteri- 
zation of cerebro vascular and tumorous diseases. Apart from the relative cerebral blood 
volume, the mean transit time is a basic diagnostic parameter. For its determination, the 
so-called impulse response function (sometimes also called impulse residue function) (21 HZj 
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of the tissue is needed, which, unfortunately, is not accessible to direct measurement. Pre- 
vious contributions consistently showed various deficiencies in regard to its determination, 
especially if diffusive tissue areas were involved. 

The aim of this contribution is to analyze such deficiencies from a more mathematical 
point of view and to point out a method that is practically feasible. Both MR and CT 
perfusion measurements are included here. 

First, we notice that the determination of the impulse response function of tissue by a 
direct deconvolution from the arterial and tissue signal does not seem to be possible in a 
reliable way. A solution is then suggested on the basis of a parametric fit. In the case of 
a disturbed blood brain barrier, it is shown that an estimation of the mean transit time 
is possible also with CT perfusion measurements, although a CT measurement is short in 
comparison to an MR measurement since CT measurement is time limited due to radiation 
exposure. 

The paper is organized as follows. First, we briefly summarize the deconvolution meth- 
ods used, and their common deficiencies in our practical situation. We continue with a 
more detailed discussion of parametric fits. Here, we discard the so-called double gamma 
fit, because it can be shown to lead to inconsistent results. Nevertheless, a closely related 
pair of parametric functions works rather well and is suggested as an alternative. 

We describe first steps towards the parametric treatment of data from defective tissue, 
and close with the suggestion of a quantity to distinguish intact from disturbed blood brain 
barrier. This is the standard deviation of the IRF and can replace the mean transit time. 
It is directly accessible even without a parametric fit, both for MR and CT data. 

2 Shortcomings of convolution techniques 

Three numerical deconvolution methods were analyzed in They are deconvolution by 

1. Fourier transform 

2. the method of moments 

3. functional calculus. 

These methods were first tested by realistic, but simulated data without noise. They 
provide acceptable and comparable results if the sampling rate is sufficiently high (which 
roughly means some five to ten sampling points per relevant interval of time). Let us 
explain the deconvolution approach in some more detail. 
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1. The deconvolution by Fourier transform uses the convolution theorem jHl p. 110], so 
that one gets 



where (.) " denotes the Fourier transform and (.) the corresponding inverse transform. 
One principle problem is the following. Whenever the Fourier transform of / and g are 
equal in a certain interval including zero, and the Fourier transform of / is zero outside 
this interval, then (/*/)"=(/* fi^) " • In general, it is thus not correct that f * gi = f * 92 
implies gi = g2- This means that the deconvolution is not well defined in general, compare 
[121 Ch. 7] for a general discussion. Fortunately, in our present context, this aspect does 
not play a major role. Another numerical problem, which is more important in our case, 
can emerge from (approximate) zeros of the denominator of ([Q). It is obvious that the 
zeros of form a subset of those of and that is always well defined from the 

theoretical point of view. However, since the Fourier transforms have to be determined 
numerically from discrete data, it can occur that the approximation of / ' features zeros 
where that of {f*g) " does not. This makes expression ((T)) undefined or leads to undesirable 
fluctuations. Finally, the numerical data for / and for f * g show similar noise level. In 
theory, f *g should be the smoother signal though, and a violation can (and actually does) 
lead to numerical sensitivity of the (discrete) inverse transform. 

2. If the function g (or a good approximation of it) is known in some parametric form, 
the unknown parameters can be determined by the method of moments. To this end, we 
need the following relation If / and g are non-negative, suitably integrable, normalized 
functions, then 



these moments exist, which is very reasonable in our context. 

Assume that it is possible to describe g hj m parameters. In this case, all moments up 
to order (m — 1) of / and f * g are needed to create a complete set of equations by means 
of Q. In our case, the functions / and f * g are non- negative, but not normalized (this 
can be mended by dividing by their respective integrals). Since, however, the moments 
have to be calculated numerically, they will show decreasing accuracy with increasing m. 
Consequently, this approach is only useful for situations with a small number of parameters. 

3. Deconvolution by the functional calculus is a less well known, but rather interesting 
method. It provides a kind of inverse formula to determine the function g if f and f * g 
are given. This does not work for any arbitrary function / but / has to be a special type 
of function. In principle, the method applies to our situation, which had been overlooked 
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in previous attempts. The following example shall illustrate this method. If 



m 




otherwise, 



t > 



(3) 



the equation (/ * g){t) = J^fit — x)g{x)dx is inverted by 



git) = {l + Dnf*gm 



(4) 



where D = is the derivation operator. This means that g can be calculated by the 
first and the second derivative of / * if / is according to 0. Further examples, and 
a systematic theory, can be found in [TT]. Unfortunately, derivatives play a major role, 
and, with increasing complexity of the convolution kernel, the derivatives required are of 
higher and higher order. Since the actual data are discrete values, the derivatives have to 
be calculated numerically. Hence, this method is suitable only to a limited extent due to 
numerical sensitivity. 

As a common problem to all three methods, the sampling rate of about one point per 
second of an MR and a CT scanner is too coarse. This is demonstrated by the following 
simple example. The function f{t) = te~*, t > 0, was convolved with a rectangular 
signal g of width 3 and height 0.2 (hence not normalized) to obtain f * g. Then, the 
particular functions were sampled at different rates. The results of deconvolution by the 
three methods are shown in Figure ^ whereas the dashed rectangle is the original function 
g for comparison. 

In comparison to this illustrative example, our real world case corresponds to the left- 
most column in Figure ^ i.e., to a clearly insufficient resolution. Note that the same 
quahtative picture would emerge for other signal types as well, not just for rectangular g. 

Additionally, the real data are superimposed by noise. From a careful inspection of our 
sample data, we estimate that noise and measurement errors together result in fiuctuations 
of five to ten percent on average. It is well known that numerical deconvolution processes 
can be very sensitive to noise. Therefore, it is preferable to describe the measured data 
by functions, e.g., by means of a parametric ansatz followed by a least squares fit. This 
requires an adequate approach. 

3 Inadequacy of the double gamma fit 

Fitting the data by gamma densities (in the literature mostly called gamma variate func- 
tions ISmni) seems to be an adequate approach. Especially the measured arterial signals 
can be described this way. But sometimes both arterial and tissue signal were displayed 
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Figure 1: Results of deconvolution by a) Fourier transform b) the method of moments and c) 
the functional calculus in dependence on the sampling rate. The latter is one point (left), two 
points (center) and ten points (right) per relevant interval of time (which is 1 in this case). For 
better illustration, the discrete sample points are successively connected by straight lines without 
additional interpolation. Note, that the function g is not normalized. 
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by gamma densities fHlEI- In this section, we want to show that such a double gamma 
fit is not reasonable. 

For simplicity, we formulate the situation in terms of probability theory. The actual 
signals are described by non-negative functions that need not be normalized. However, the 
normalization constants play no role here (one can simply divide by them to get normalized 
signals), so that we can profit from the standard tools of probability theory. 

If / is an integrable non-negative function. 



(f)[t) = / e''''f{x)dx (5) 

is the (inverse) Fourier transform of /. If / is a probability density, is called its charac- 
teristic function, see jHj for details. In the case of / = fx^a being a (normalized) gamma 
density, i.e., 

, , fo, X < 0, 

oo 

with A > 0, a > 1 and T{a) = J x"^~^e~^dx, the characteristic function is [HI p. 343] 



0A,a(t) = ' ^^^^ ^ ^ ^' 

which can be calculated with elementary methods from complex analysis. 

Starting from two gamma densities, we consider the convolution /a^q * 9 = f'y,u- (In 
doing so, fx^a corresponds to the measured arterial signal and f^^^, to the tissue signal, 
whereas g is the IRF wanted, up to a multiplicative factor). By means of the convolution 
theorem and the characteristic function of the gamma density, one gets the characteristic 
function of g as 

0^(t) = f =: h{t). (8) 

The central question now is whether h is the characteristic function of a probability density. 
An answer is given by the Bochner-Khintchine theorem which states that a continuous 
function ip{t), t G M, with tf{0) = 1 is the characteristic function of a probability measure 
if and only if it is positive definite, see the appendix for details. 

The violation of some elementary properties of positive definite functions in this case 
already implies that h is not positive definite in general. Fitting the data by gamma 
densities provides parameters in a region where property (3) of the appendix is violated, 
see for details. Consequently, h is not a characteristic function and therefore g is not 
a probability density. This, in turn, means that g takes negative values which cannot be 
interpreted in our context. Consequently, the sometimes used double gamma fit is not 
consistent in this context and should be avoided. 
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4 An alternative approach 



Let us look for an alternative. It is reasonable (and also rather plausible) to assume the 
IRF itself to be gamma distributed. In this case, we can calculate the shape of the tissue 
signal if the arterial signal is still described by a gamma density. We have to determine 
f * g provided that / = fx^a and g = g^^p according to Once again, by means of the 
convolution theorem and the characteristic function of the gamma density, one now finds 

with 0/=i,g(O) = 1. The corresponding inverse transform reads jTOit Ch. 3.2] 

fs\ (10) 

= — - — — — iFi[a; a + p;{fi- X)x] 
L[a + p) 

for X > and {f * g){x) = otherwise, where iFi is the confluent hypergeometric function 
P Ch. 13] with 

whereas (c)^ := c(c + 1) • . . . ■ (c + A; — 1) is the Pochhammer symbol. The normalization 
of the inverse transform is chosen such that (/ * g){x) is indeed a normalized probability 
density, i.e., (/ * g){x) > together with (/ * g){x)dx = 1. 

The two formulations of (|T0|) are equal by Kummer's functional equation for iFi, see 
[H Eq. 13.1.27]. For A > /i (resp. A < /i) the first (resp. the second) version immediately 
shows the non- negativity of (/ * g){x). For X = fi, f * g is the density of a gamma 
distribution, in hne with the convolution semigroup property (/a,q * f\,(3){x) = fx^a+p^x) 
for gamma densities with the same parameter A, compare 0. 

The question is now whether the tissue signal can be described by the expression in (^UJ. 
An example (compare with Figure Ej) shows that the fit is rather convincing, so that this 
ansatz is phenomenologically sound. Of course, it remains an interesting open question at 
this point whether further corroboration, e.g., by microscopic modelling, is possible. 

The function in Eq. (|1U|) fits the tissue signal at least as well as a gamma density would 
do. But in contrast to the latter, we no longer encounter problems with negative values of 
the IRF. 

Note, however, that this approach can only be used for tissue with an intact blood 
brain barrier because the fit separates the recirculation from the real signal. Tissue with 
a defective blood brain barrier, such as a brain tumor, has to be considered separately 
because recirculation is superimposed by diffusion. 
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Figure 2: Parametric fit of tlie tissue signal measured by an MR scanner (top) and a CT scanner 
(bottom), witli tlie function from Eq. H1U|) . up to an additional normalization factor. All points 
to the left of the dashed vertical line were used for the fit. 



4.1 Tissue with defective blood brain barrier 

The measured signal of a tissue with a permeable blood brain barrier (e.g., brain tumor) 
consists of a diffusive and a perfusive part |3IH]- The idea here is to split up these parts and 
deconvolve them separately. This way, additional features of the different contributions 
can be detected and displayed more clearly. 

After injection of the contrast medium, most of these molecules remain in the intravas- 
cular space during their transit through the vascular network. However, a small portion 
diffuses into the extravascular space and then back into the intravascular space due to the 
defective blood brain barrier. Because diffusion is much slower than perfusion, contrast 
medium molecules still diffuse into the extravascular space during recirculation El HE] • 

To our knowledge, there is no theoretical derivation known to describe the diffusive 
part in detail, wherefore a phenomenological approach is necessary and appropriate. The 
function 



with a>0, 6>l,c>0 and d > unknown parameters, seems to be applicable to capture 
the diffusion part parametrically, see Figure El 

The perfusive part can be obtained approximately by subtraction of the diffusion part 
from the original signal (compare with Figure E] top). Figure 0] bottom shows that (|1(J|) 
leads to a fit of acceptable quality. This means that the method introduced above can also 
be applied to determine the IRF of the perfusion. In this case, the IRF is a gamma density 
with the associated parameters from the fit. 

It remains to consider the diffusive part. We did not find an explicit and practically use- 
ful expression for the IRF, so that we employed a numerical deconvolution technique (e.g., 
by Fourier transform). This, in contrast to above, does not pose a problem because both 
signals required are now available in parametrized form, and because the non-uniqueness 
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Figure 3: The function f{t) = at^e '^^'^ as a parametric fit for the diffusive part. All points to 
the right of the dashed vertical line were used to determine the associated parameters. 
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Figure 4: Perfusion after subtraction of the diffusive part (top) and Eq. l\W\i as an adequate fit 
of the perfusion (bottom). All points to the left of the dashed line were used in this fit. The few 
points to the right not fitted by (|lfl|) are the remainder of recirculation. 



problem does not show up. Now, one can put the IRFs together to obtain the IRF of 
tissue with a permeable blood brain barrier. Figure El shows the IRFs for the different 
tissue types. On the top, the IRF for tissue with an intact blood brain barrier is shown, 
and the IRF for tissue with a defective blood brain barrier on the bottom. 
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Figure 5: The IRF for tissue with an intact blood brain barrier (top) and a defective blood brain 
barrier (bottom). Note the different scale on the time axis for comparison. 
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4.2 Diagnostic parameters 

So far, the area to height relation has been used to determine the mean transit time [MTT) 
[3 EI- There, the area under the curve of the IRF is divided by the maximum of the IRF, 
i.e., 

f R(t)dt , , 

MTT = ^ — (12) 

Rmax 

if the IRF is denoted by R. However, this approach is only apphcable for tissue with an 
intact blood brain barrier under the simplifying (and somewhat unrealistic) assumption 
that the IRF is a pulse of rectangular shape. By application of this formula to tissue with 
a permeable blood brain barrier, a rather unrealistic and unreliable estimate of the mean 
transit time is obtained. For characterizing and possibly detecting a brain tumor, a more 
reliable quantity is required. 

The standard deviation of the IRF appears to be an adequate measure for the mean 
transit time. To expand on this suggestion, let / be a non-negative signal with a compact 
support. If = J^f{t)dt is the normalization constant. 



ruf 



^him (13) 



is the signal mean and 



is its standard deviation. 

The standard deviation features the correct dimension, namely time, and does not 
change if the entire signal is shifted along the time axis. More importantly, its definition 
is independent of the shape of the IRF in the sense that it does not require any specific 
assumptions on the signal form. Consequently, it is applicable to the different tissue types. 

In the case of the IRF being a non-normalized gamma density, there is an asymptotically 
linear relation between the mean transit time, computed by the area to height formula, 
and the standard deviation of the IRF. Thus, it is no big deal to replace the traditional 
formula by the standard deviation, because one calculates nothing much different in this 
case. 

The theoretical relation can be derived as follows. Let the IRF again be denoted by R 
and assume that it is a non-normalized gamma density with parameters b, (3 and /x, i.e., 

R{t) = { - (15) 

0, otherwise. 
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The standard deviation reads ar = the area under the curve is J^R{t)dt = b and 

it is easily checked that the maximum value is Rmax = = ^r^{^^)^ ^e^"^. 

Calculating the ratio of err and MTTr from (fT^ leads to 

err _ VPiP- ir^' (^g) 



MTTr r(/3) e^-i ' 

This expression depends on parameter j3 only. Figure IHl shows this function in dependence 
on p. The dashed horizontal line indicates the asymptotic value of the function as — > cxd. 



a/MTT 

Ir 



Figure 6: Theoretical relation between a and MTT in the case of a gamma density, in dependence 
on the parameter (3. The dashed line runs at height ~ 0.39. 

V 27r 

Stirling's formula [TJ Eq. 6.1.37] states that 

T{x) = e"" 72^ (1 + Oix'')) (17) 
as X — > oo. With this expression we get 

^ ^ +0{(3-') ^ ^^0.39. (18) 



MTTr (3-1 " ' ^ 

Using the standard deviation of the IRF as a measure for the mean transit time, one gets 
significantly different values in dependence on the tissue type. Tissue with a permeable 
blood brain barrier (e.g., brain tumor) possesses a greater mean transit time than normal 
tissue (compare with figure El again) . 

To turn this into a practical method, the following ansatz looks promising. It is possible 
to determine the standard deviation of the IRF directly from the given (noisy) data without 
making an explicit deconvolution. Equation (j2I) implies ctj^^ = aj + ag, if / and g are 
probability densities. Consequently, the standard deviation of g can be calculated as 



^. = \R.-^/- (19) 
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In our case, we have to normalize the signals prior to applying ()19|). This can be done 
on the basis of the discrete data directly. Our signals are available as a finite time series 
{he I 1 < < Timax}-, where we assume, for simplicity, that the integral index i is the time 
(this being literally true for our MR test data). The normalization leads to 



where are the non-normalized values of the respective signal and pk the normalized ones. 
This way, one gets discrete probability distributions. Therefore, mean m and variance cr^ 
read 



Using Eq. (fTI?|) . the standard deviation of the IRF can now be computed. 

Taking the values from the first pass of the signals (which refers to their structure in 
time, see jS]), the standard deviation provides already differences in dependence on the 
tissue type. We tested this with some exemplary patient data. For both MR and CT data, 
the values were between 1.8 and 2.6 in the case of normal tissue and between 3.1 and 3.4 
in the case of a brain tumor. This way, getting a first indication in respect of the tissue 
type seems possible. Amazingly, such information can also be extracted from CT data, in 
spite of the rather short time series available. 

5 Concluding remarks 

Our brief exposition shows two things. First, a reasonable parametric description of blood 
fiow data, both from MR and CT, is possible and does not suffer from standard problems 
of numerical deconvolution methods. 

Second, the IRF is sufficiently accessible to draw first conclusions on the tissue type 
from it. In this context, we suggest the use of the standard deviation of the IRF as a 
potential diagnostic parameter. 

All our methods were tested on real data, which are contained in However, this 
clearly needs further investigation, and a systematic test with a much larger data set. We 
hope to report on some progress in the near future. 

Appendix: Positivity and positive definiteness 

Let us explain the connection between these concepts, where we follow A complex- 
valued function ip on the real line is called positive definite, if for all n G N and for all 




(21) 
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n-tupels {xi, . . . , Xn) with Xi G R, the nxn matrix ( (f{xi — Xj) )i<i,j<n is positive Hermitian. 
This definition imphes the following properties: 



(1) cpiO) > 

(2) (p{-x) = Lp{x), for all X G M 

(3) \^{x)\<^{0), forallxGM 

(4) \ifix)~^{y)\'< 2^{0) {^{0)-Re{ip{x-y))), 

see jU Ch. 3.4] for details. Property (4) is often called Krein's inequality. 

A function tp is positive definite if and only if there is a finite positive measure on M 



Moreover, is real-valued if the measure /i is symmetric, i.e., fi{A) = fi{—A) for all Borel 
subsets of R. This is a special case of Bochner's theorem, compare [H Thm. 3.12], and fi 
is uniquely determined. 

The special case that is a probability measure leads to ip{0) = 1. The function is now 
the characteristic function of n, and the correspondence between positive definiteness of 
if, with (p{0) = 1, and /z being a probability measure is also called the Bochner-Khintchine 
theorem jl4j . 

In our case, the function ip is given and we need to know whether the corresponding /i, 
the IRF, is a positive signal function. The latter is a special case of a positive measure. 

If, however, any of the properties (1) - (4) fails, we can conclude that cannot be a 
positive measure, and hence certainly not a positive signal. 
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